home *** CD-ROM | disk | FTP | other *** search
/ EuroCD 3 / EuroCD 3.iso / Programming / vbcc / machines / amiga68k / libsrc / math / math_040 / log.s < prev    next >
Encoding:
Text File  |  1998-06-24  |  24.5 KB  |  522 lines
  1. *
  2. *   $VER: log.s 33.1 (20.1.97)
  3. *
  4. *   " computes the natural logarithm of a normalized input"
  5. *
  6. *   Version history:
  7. *
  8. *   33.1    20.1.97 (c) Motorola
  9. *
  10. *           - cut'n pasted from M68060SP
  11. *
  12. *   33.2    22.1.97 (c) Motorola
  13. *
  14. *           - added log10 and log2
  15. *           - log trashed a6, fixed
  16. *
  17.  
  18.         machine 68040
  19.         fpu     1
  20.  
  21.         XDEF    _log
  22.         XDEF    @log
  23.         XDEF    _log10
  24.         XDEF    @log10
  25.         XDEF    _log2
  26.         XDEF    @log2
  27.  
  28. *************************************************************************
  29. * log():    computes the natural logarithm of a normalized input        *
  30. *                                                                       *
  31. * INPUT *************************************************************** *
  32. *       fp0 = extended precision input                                  *
  33. *                                                                       *
  34. * OUTPUT ************************************************************** *
  35. *       fp0 = log(X)                                                    *
  36. *                                                                       *
  37. * ACCURACY and MONOTONICITY ******************************************* *
  38. *       The returned result is within 2 ulps in 64 significant bit,     *
  39. *       i.e. within 0.5001 ulp to 53 bits if the result is subsequently *
  40. *       rounded to double precision. The result is provably monotonic   *
  41. *       in double precision.                                            *
  42. *                                                                       *
  43. * ALGORITHM *********************************************************** *
  44. *       LOGN:                                                           *
  45. *       Step 1. If |X-1| < 1/16, approximate log(X) by an odd           *
  46. *               polynomial in u, where u = 2(X-1)/(X+1). Otherwise,     *
  47. *               move on to Step 2.                                      *
  48. *                                                                       *
  49. *       Step 2. X = 2**k * Y where 1 <= Y < 2. Define F to be the first *
  50. *               seven significant bits of Y plus 2**(-7), i.e.          *
  51. *               F = 1.xxxxxx1 in base 2 where the six "x" match those   *
  52. *               of Y. Note that |Y-F| <= 2**(-7).                       *
  53. *                                                                       *
  54. *       Step 3. Define u = (Y-F)/F. Approximate log(1+u) by a           *
  55. *               polynomial in u, log(1+u) = poly.                       *
  56. *                                                                       *
  57. *       Step 4. Reconstruct                                             *
  58. *               log(X) = log( 2**k * Y ) = k*log(2) + log(F) + log(1+u) *
  59. *               by k*log(2) + (log(F) + poly). The values of log(F) are *
  60. *               calculated beforehand and stored in the program.        *
  61. *                                                                       *
  62. *       Implementation Notes:                                           *
  63. *       Note 1. There are 64 different possible values for F, thus 64   *
  64. *               log(F)'s need to be tabulated. Moreover, the values of  *
  65. *               1/F are also tabulated so that the division in (Y-F)/F  *
  66. *               can be performed by a multiplication.                   *
  67. *                                                                       *
  68. *       Note 2. To fully exploit the pipeline, polynomials are usually  *
  69. *               separated into two parts evaluated independently before *
  70. *               being added up.                                         *
  71. *                                                                       *
  72. *************************************************************************
  73.  
  74.         cnop            0,4
  75.  
  76. LOGOF2  dc.l            $3FFE0000,$B17217F7,$D1CF79AC,$00000000
  77. one     dc.l            $3F800000
  78. zero    dc.l            $00000000
  79. infty   dc.l            $7F800000
  80. negone  dc.l            $BF800000
  81. LOGA6   dc.l            $3FC2499A,$B5E4040B
  82. LOGA5   dc.l            $BFC555B5,$848CB7DB
  83. LOGA4   dc.l            $3FC99999,$987D8730
  84. LOGA3   dc.l            $BFCFFFFF,$FF6F7E97
  85. LOGA2   dc.l            $3FD55555,$555555A4
  86. LOGA1   dc.l            $BFE00000,$00000008
  87. LOGB5   dc.l            $3F175496,$ADD7DAD6
  88. LOGB4   dc.l            $3F3C71C2,$FE80C7E0
  89. LOGB3   dc.l            $3F624924,$928BCCFF
  90. LOGB2   dc.l            $3F899999,$999995EC
  91. LOGB1   dc.l            $3FB55555,$55555555
  92. TWO     dc.l            $40000000,$00000000
  93. LTHOLD  dc.l            $3f990000,$80000000,$00000000,$00000000
  94. LOGTBL  dc.l            $3FFE0000,$FE03F80F,$E03F80FE,$00000000
  95.         dc.l            $3FF70000,$FF015358,$833C47E2,$00000000
  96.         dc.l            $3FFE0000,$FA232CF2,$52138AC0,$00000000
  97.         dc.l            $3FF90000,$BDC8D83E,$AD88D549,$00000000
  98.         dc.l            $3FFE0000,$F6603D98,$0F6603DA,$00000000
  99.         dc.l            $3FFA0000,$9CF43DCF,$F5EAFD48,$00000000
  100.         dc.l            $3FFE0000,$F2B9D648,$0F2B9D65,$00000000
  101.         dc.l            $3FFA0000,$DA16EB88,$CB8DF614,$00000000
  102.         dc.l            $3FFE0000,$EF2EB71F,$C4345238,$00000000
  103.         dc.l            $3FFB0000,$8B29B775,$1BD70743,$00000000
  104.         dc.l            $3FFE0000,$EBBDB2A5,$C1619C8C,$00000000
  105.         dc.l            $3FFB0000,$A8D839F8,$30C1FB49,$00000000
  106.         dc.l            $3FFE0000,$E865AC7B,$7603A197,$00000000
  107.         dc.l            $3FFB0000,$C61A2EB1,$8CD907AD,$00000000
  108.         dc.l            $3FFE0000,$E525982A,$F70C880E,$00000000
  109.         dc.l            $3FFB0000,$E2F2A47A,$DE3A18AF,$00000000
  110.         dc.l            $3FFE0000,$E1FC780E,$1FC780E2,$00000000
  111.         dc.l            $3FFB0000,$FF64898E,$DF55D551,$00000000
  112.         dc.l            $3FFE0000,$DEE95C4C,$A037BA57,$00000000
  113.         dc.l            $3FFC0000,$8DB956A9,$7B3D0148,$00000000
  114.         dc.l            $3FFE0000,$DBEB61EE,$D19C5958,$00000000
  115.         dc.l            $3FFC0000,$9B8FE100,$F47BA1DE,$00000000
  116.         dc.l            $3FFE0000,$D901B203,$6406C80E,$00000000
  117.         dc.l            $3FFC0000,$A9372F1D,$0DA1BD17,$00000000
  118.         dc.l            $3FFE0000,$D62B80D6,$2B80D62C,$00000000
  119.         dc.l            $3FFC0000,$B6B07F38,$CE90E46B,$00000000
  120.         dc.l            $3FFE0000,$D3680D36,$80D3680D,$00000000
  121.         dc.l            $3FFC0000,$C3FD0329,$06488481,$00000000
  122.         dc.l            $3FFE0000,$D0B69FCB,$D2580D0B,$00000000
  123.         dc.l            $3FFC0000,$D11DE0FF,$15AB18CA,$00000000
  124.         dc.l            $3FFE0000,$CE168A77,$25080CE1,$00000000
  125.         dc.l            $3FFC0000,$DE1433A1,$6C66B150,$00000000
  126.         dc.l            $3FFE0000,$CB8727C0,$65C393E0,$00000000
  127.         dc.l            $3FFC0000,$EAE10B5A,$7DDC8ADD,$00000000
  128.         dc.l            $3FFE0000,$C907DA4E,$871146AD,$00000000
  129.         dc.l            $3FFC0000,$F7856E5E,$E2C9B291,$00000000
  130.         dc.l            $3FFE0000,$C6980C69,$80C6980C,$00000000
  131.         dc.l            $3FFD0000,$82012CA5,$A68206D7,$00000000
  132.         dc.l            $3FFE0000,$C4372F85,$5D824CA6,$00000000
  133.         dc.l            $3FFD0000,$882C5FCD,$7256A8C5,$00000000
  134.         dc.l            $3FFE0000,$C1E4BBD5,$95F6E947,$00000000
  135.         dc.l            $3FFD0000,$8E44C60B,$4CCFD7DE,$00000000
  136.         dc.l            $3FFE0000,$BFA02FE8,$0BFA02FF,$00000000
  137.         dc.l            $3FFD0000,$944AD09E,$F4351AF6,$00000000
  138.         dc.l            $3FFE0000,$BD691047,$07661AA3,$00000000
  139.         dc.l            $3FFD0000,$9A3EECD4,$C3EAA6B2,$00000000
  140.         dc.l            $3FFE0000,$BB3EE721,$A54D880C,$00000000
  141.         dc.l            $3FFD0000,$A0218434,$353F1DE8,$00000000
  142.         dc.l            $3FFE0000,$B92143FA,$36F5E02E,$00000000
  143.         dc.l            $3FFD0000,$A5F2FCAB,$BBC506DA,$00000000
  144.         dc.l            $3FFE0000,$B70FBB5A,$19BE3659,$00000000
  145.         dc.l            $3FFD0000,$ABB3B8BA,$2AD362A5,$00000000
  146.         dc.l            $3FFE0000,$B509E68A,$9B94821F,$00000000
  147.         dc.l            $3FFD0000,$B1641795,$CE3CA97B,$00000000
  148.         dc.l            $3FFE0000,$B30F6352,$8917C80B,$00000000
  149.         dc.l            $3FFD0000,$B7047551,$5D0F1C61,$00000000
  150.         dc.l            $3FFE0000,$B11FD3B8,$0B11FD3C,$00000000
  151.         dc.l            $3FFD0000,$BC952AFE,$EA3D13E1,$00000000
  152.         dc.l            $3FFE0000,$AF3ADDC6,$80AF3ADE,$00000000
  153.         dc.l            $3FFD0000,$C2168ED0,$F458BA4A,$00000000
  154.         dc.l            $3FFE0000,$AD602B58,$0AD602B6,$00000000
  155.         dc.l            $3FFD0000,$C788F439,$B3163BF1,$00000000
  156.         dc.l            $3FFE0000,$AB8F69E2,$8359CD11,$00000000
  157.         dc.l            $3FFD0000,$CCECAC08,$BF04565D,$00000000
  158.         dc.l            $3FFE0000,$A9C84A47,$A07F5638,$00000000
  159.         dc.l            $3FFD0000,$D2420487,$2DD85160,$00000000
  160.         dc.l            $3FFE0000,$A80A80A8,$0A80A80B,$00000000
  161.         dc.l            $3FFD0000,$D7894992,$3BC3588A,$00000000
  162.         dc.l            $3FFE0000,$A655C439,$2D7B73A8,$00000000
  163.         dc.l            $3FFD0000,$DCC2C4B4,$9887DACC,$00000000
  164.         dc.l            $3FFE0000,$A4A9CF1D,$96833751,$00000000
  165.         dc.l            $3FFD0000,$E1EEBD3E,$6D6A6B9E,$00000000
  166.         dc.l            $3FFE0000,$A3065E3F,$AE7CD0E0,$00000000
  167.         dc.l            $3FFD0000,$E70D785C,$2F9F5BDC,$00000000
  168.         dc.l            $3FFE0000,$A16B312E,$A8FC377D,$00000000
  169.         dc.l            $3FFD0000,$EC1F392C,$5179F283,$00000000
  170.         dc.l            $3FFE0000,$9FD809FD,$809FD80A,$00000000
  171.         dc.l            $3FFD0000,$F12440D3,$E36130E6,$00000000
  172.         dc.l            $3FFE0000,$9E4CAD23,$DD5F3A20,$00000000
  173.         dc.l            $3FFD0000,$F61CCE92,$346600BB,$00000000
  174.         dc.l            $3FFE0000,$9CC8E160,$C3FB19B9,$00000000
  175.         dc.l            $3FFD0000,$FB091FD3,$8145630A,$00000000
  176.         dc.l            $3FFE0000,$9B4C6F9E,$F03A3CAA,$00000000
  177.         dc.l            $3FFD0000,$FFE97042,$BFA4C2AD,$00000000
  178.         dc.l            $3FFE0000,$99D722DA,$BDE58F06,$00000000
  179.         dc.l            $3FFE0000,$825EFCED,$49369330,$00000000
  180.         dc.l            $3FFE0000,$9868C809,$868C8098,$00000000
  181.         dc.l            $3FFE0000,$84C37A7A,$B9A905C9,$00000000
  182.         dc.l            $3FFE0000,$97012E02,$5C04B809,$00000000
  183.         dc.l            $3FFE0000,$87224C2E,$8E645FB7,$00000000
  184.         dc.l            $3FFE0000,$95A02568,$095A0257,$00000000
  185.         dc.l            $3FFE0000,$897B8CAC,$9F7DE298,$00000000
  186.         dc.l            $3FFE0000,$94458094,$45809446,$00000000
  187.         dc.l            $3FFE0000,$8BCF55DE,$C4CD05FE,$00000000
  188.         dc.l            $3FFE0000,$92F11384,$0497889C,$00000000
  189.         dc.l            $3FFE0000,$8E1DC0FB,$89E125E5,$00000000
  190.         dc.l            $3FFE0000,$91A2B3C4,$D5E6F809,$00000000
  191.         dc.l            $3FFE0000,$9066E68C,$955B6C9B,$00000000
  192.         dc.l            $3FFE0000,$905A3863,$3E06C43B,$00000000
  193.         dc.l            $3FFE0000,$92AADE74,$C7BE59E0,$00000000
  194.         dc.l            $3FFE0000,$8F1779D9,$FDC3A219,$00000000
  195.         dc.l            $3FFE0000,$94E9BFF6,$15845643,$00000000
  196.         dc.l            $3FFE0000,$8DDA5202,$37694809,$00000000
  197.         dc.l            $3FFE0000,$9723A1B7,$20134203,$00000000
  198.         dc.l            $3FFE0000,$8CA29C04,$6514E023,$00000000
  199.         dc.l            $3FFE0000,$995899C8,$90EB8990,$00000000
  200.         dc.l            $3FFE0000,$8B70344A,$139BC75A,$00000000
  201.         dc.l            $3FFE0000,$9B88BDAA,$3A3DAE2F,$00000000
  202.         dc.l            $3FFE0000,$8A42F870,$5669DB46,$00000000
  203.         dc.l            $3FFE0000,$9DB4224F,$FFE1157C,$00000000
  204.         dc.l            $3FFE0000,$891AC73A,$E9819B50,$00000000
  205.         dc.l            $3FFE0000,$9FDADC26,$8B7A12DA,$00000000
  206.         dc.l            $3FFE0000,$87F78087,$F78087F8,$00000000
  207.         dc.l            $3FFE0000,$A1FCFF17,$CE733BD4,$00000000
  208.         dc.l            $3FFE0000,$86D90544,$7A34ACC6,$00000000
  209.         dc.l            $3FFE0000,$A41A9E8F,$5446FB9F,$00000000
  210.         dc.l            $3FFE0000,$85BF3761,$2CEE3C9B,$00000000
  211.         dc.l            $3FFE0000,$A633CD7E,$6771CD8B,$00000000
  212.         dc.l            $3FFE0000,$84A9F9C8,$084A9F9D,$00000000
  213.         dc.l            $3FFE0000,$A8489E60,$0B435A5E,$00000000
  214.         dc.l            $3FFE0000,$83993052,$3FBE3368,$00000000
  215.         dc.l            $3FFE0000,$AA59233C,$CCA4BD49,$00000000
  216.         dc.l            $3FFE0000,$828CBFBE,$B9A020A3,$00000000
  217.         dc.l            $3FFE0000,$AC656DAE,$6BCC4985,$00000000
  218.         dc.l            $3FFE0000,$81848DA8,$FAF0D277,$00000000
  219.         dc.l            $3FFE0000,$AE6D8EE3,$60BB2468,$00000000
  220.         dc.l            $3FFE0000,$80808080,$80808081,$00000000
  221.         dc.l            $3FFE0000,$B07197A2,$3C46C654,$00000000
  222.  
  223.  
  224. X           EQU         -12
  225. XDCARE      EQU         X+2
  226. XFRAC       EQU         X+4
  227. KLOG2       EQU         -12
  228. SAVEU       EQU         -12
  229. F           EQU         -24
  230. FFRAC       EQU         F+4
  231. TEMP_SIZE   EQU         24
  232.  
  233. _log
  234.         fmove.d         (4,sp),fp0
  235. @log
  236.         link            a1,#-TEMP_SIZE
  237.         fmove.x         fp0,(X,a1)
  238.  
  239. .LOGBGN
  240. ;--FPCR SAVED AND CLEARED, INPUT IS 2^(ADJK)*FP0, FP0 CONTAINS
  241. ;--A FINITE, NON-ZERO, NORMALIZED NUMBER.
  242.  
  243.         move.l          (X,a1),d1
  244.         move.w          (X+4,a1),d1
  245.  
  246.         tst.l           d1                      ; CHECK IF X IS NEGATIVE
  247.         ble.w           .LOGNEG                 ; LOG OF NEGATIVE ARGUMENT IS INVALID
  248.  
  249. ; X IS POSITIVE, CHECK IF X IS NEAR 1
  250.  
  251.         cmp.l           #$3ffef07d,d1           ; IS X < 15/16?
  252.         blt.b           .LOGMAIN                ; YES
  253.         cmp.l           #$3fff8841,d1           ; IS X > 17/16?
  254.         ble.w           .LOGNEAR1               ; NO
  255.  
  256. .LOGMAIN
  257.  
  258. ;--THIS SHOULD BE THE USUAL CASE, X NOT VERY CLOSE TO 1
  259.  
  260. ;--X = 2^(K) * Y, 1 <= Y < 2. THUS, Y = 1.XXXXXXXX....XX IN BINARY.
  261. ;--WE DEFINE F = 1.XXXXXX1, I.E. FIRST 7 BITS OF Y AND ATTACH A 1.
  262. ;--THE IDEA IS THAT LOG(X) = K*LOG2 + LOG(Y)
  263. ;--                      = K*LOG2 + LOG(F) + LOG(1 + (Y-F)/F).
  264. ;--NOTE THAT U = (Y-F)/F IS VERY SMALL AND THUS APPROXIMATING
  265. ;--LOG(1+U) CAN BE VERY EFFICIENT.
  266. ;--ALSO NOTE THAT THE VALUE 1/F IS STORED IN A TABLE SO THAT NO
  267. ;--DIVISION IS NEEDED TO CALCULATE (Y-F)/F.
  268.  
  269. ;--GET K, Y, F, AND ADDRESS OF 1/F.
  270.  
  271.         asr.l           #8,d1
  272.         asr.l           #8,d1                   ; SHIFTED 16 BITS, BIASED EXPO. OF X
  273.         sub.l           #$3FFF,d1               ; THIS IS K
  274.  
  275.         lea             (LOGTBL,pc),a0          ; BASE ADDRESS OF 1/F AND LOG(F)
  276.         fmove.l         d1,fp1                  ; CONVERT K TO FLOATING-POINT FORMAT
  277.  
  278. ;--WHILE THE CONVERSION IS GOING ON, WE GET F AND ADDRESS OF 1/F
  279.         move.l          #$3FFF0000,(X,a1)       ; X IS NOW Y, I.E. 2^(-K)*X
  280.         move.l          (XFRAC,a1),(FFRAC,a1)
  281.         and.l           #$FE000000,(FFRAC,a1)   ; FIRST 7 BITS OF Y
  282.         or.l            #$01000000,(FFRAC,a1)
  283.         move.l          (FFRAC,a1),d1
  284.         and.l           #$7E000000,d1
  285.         asr.l           #8,d1
  286.         asr.l           #8,d1
  287.         asr.l           #4,d1                   ; SHIFTED 20, D0 IS THE DISPLACEMENT
  288.         add.l           d1,a0                   ; A0 IS THE ADDRESS FOR 1/F
  289.  
  290.         fmove.x         (X,a1),fp0
  291.         move.l          #$3fff0000,(F,a1)
  292.         clr.l           (F+8,a1)
  293.         fsub.x          (F,a1),fp0              ; Y-F
  294.         fmovem.x        fp2/fp3,-(sp)           ; SAVE FP2-3 WHILE FP0 IS NOT READY
  295. ;--SUMMARY: FP0 IS Y-F, A0 IS ADDRESS OF 1/F, FP1 IS K
  296. ;--REGISTERS SAVED: FPCR, FP1, FP2
  297.  
  298. .LP1CONT1
  299. ;--AN RE-ENTRY POINT FOR LOGNP1
  300.         fmul.x          (a0),fp0                ; FP0 IS U = (Y-F)/F
  301.         fmul.x          (LOGOF2,pc),fp1         ; GET K*LOG2 WHILE FP0 IS NOT READY
  302.         fmove.x         fp0,fp2
  303.         fmul.x          fp2,fp2                 ; FP2 IS V=U*U
  304.         fmove.x         fp1,(KLOG2,a1)          ; PUT K*LOG2 IN MEMEORY, FREE FP1
  305.  
  306. ;--LOG(1+U) IS APPROXIMATED BY
  307. ;--U + V*(A1+U*(A2+U*(A3+U*(A4+U*(A5+U*A6))))) WHICH IS
  308. ;--[U + V*(A1+V*(A3+V*A5))]  +  [U*V*(A2+V*(A4+V*A6))]
  309.  
  310.         fmove.x         fp2,fp3
  311.         fmove.x         fp2,fp1
  312.  
  313.         fmul.d          (LOGA6,pc),fp1          ; V*A6
  314.         fmul.d          (LOGA5,pc),fp2          ; V*A5
  315.  
  316.         fadd.d          (LOGA4,pc),fp1          ; A4+V*A6
  317.         fadd.d          (LOGA3,pc),fp2          ; A3+V*A5
  318.  
  319.         fmul.x          fp3,fp1                 ; V*(A4+V*A6)
  320.         fmul.x          fp3,fp2                 ; V*(A3+V*A5)
  321.  
  322.         fadd.d          (LOGA2,pc),fp1          ; A2+V*(A4+V*A6)
  323.         fadd.d          (LOGA1,pc),fp2          ; A1+V*(A3+V*A5)
  324.  
  325.         fmul.x          fp3,fp1                 ; V*(A2+V*(A4+V*A6))
  326.         addq.l          #8,a0                   ; ADDRESS OF LOG(F)
  327.         fmul.x          fp3,fp2                 ; V*(A1+V*(A3+V*A5))
  328.         addq.l          #8,a0
  329.         fmul.x          fp0,fp1                 ; U*V*(A2+V*(A4+V*A6))
  330.         fadd.x          fp2,fp0                 ; U+V*(A1+V*(A3+V*A5))
  331.  
  332.         fadd.x          (a0),fp1                ; LOG(F)+U*V*(A2+V*(A4+V*A6))
  333.         fmovem.x        (sp)+,fp2/fp3           ; RESTORE FP2-3
  334.         fadd.x          fp1,fp0                 ; FP0 IS LOG(F) + LOG(1+U)
  335.  
  336.         fadd.x          (KLOG2,a1),fp0          ; FINAL ADD
  337.         unlk            a1
  338.         rts
  339.  
  340. .LOGNEAR1
  341.  
  342. ; if the input is exactly equal to one, then exit through ld_pzero.
  343. ; if these 2 lines weren't here, the correct answer would be returned
  344. ; but the INEX2 bit would be set.
  345.  
  346.         fcmp.s          #1,fp0                  ; is it equal to one?
  347.         fbeq            .pzero                  ; yes
  348.  
  349. ;--REGISTERS SAVED: FPCR, FP1. FP0 CONTAINS THE INPUT.
  350.         fmove.x         fp0,fp1
  351.         fsub.s          (one,pc),fp1            ; FP1 IS X-1
  352.         fadd.s          (one,pc),fp0            ; FP0 IS X+1
  353.         fadd.x          fp1,fp1                 ; FP1 IS 2(X-1)
  354. ;--LOG(X) = LOG(1+U/2)-LOG(1-U/2) WHICH IS AN ODD POLYNOMIAL
  355. ;--IN U, U = 2(X-1)/(X+1) = FP1/FP0
  356.  
  357. .LP1CONT2
  358. ;--THIS IS AN RE-ENTRY POINT FOR LOGNP1
  359.         fdiv.x          fp0,fp1                 ; FP1 IS U
  360.         fmovem.x        fp2/fp3,-(sp)           ; SAVE FP2-3
  361.  
  362. ;--REGISTERS SAVED ARE NOW FPCR,FP1,FP2,FP3
  363. ;--LET V=U*U, W=V*V, CALCULATE
  364. ;--U + U*V*(B1 + V*(B2 + V*(B3 + V*(B4 + V*B5)))) BY
  365. ;--U + U*V*(  [B1 + W*(B3 + W*B5)]  +  [V*(B2 + W*B4)]  )
  366.  
  367.         fmove.x         fp1,fp0
  368.         fmul.x          fp0,fp0                 ; FP0 IS V
  369.         fmove.x         fp1,(SAVEU,a1)          ; STORE U IN MEMORY, FREE FP1
  370.         fmove.x         fp0,fp1
  371.         fmul.x          fp1,fp1                 ; FP1 IS W
  372.  
  373.         fmove.d         (LOGB5,pc),fp3
  374.         fmove.d         (LOGB4,pc),fp2
  375.  
  376.         fmul.x          fp1,fp3                 ; W*B5
  377.         fmul.x          fp1,fp2                 ; W*B4
  378.  
  379.         fadd.d          (LOGB3,pc),fp3          ; B3+W*B5
  380.         fadd.d          (LOGB2,pc),fp2          ; B2+W*B4
  381.  
  382.         fmul.x          fp3,fp1                 ; W*(B3+W*B5), FP3 RELEASED
  383.  
  384.         fmul.x          fp0,fp2                 ; V*(B2+W*B4)
  385.  
  386.         fadd.d          (LOGB1,pc),fp1          ; B1+W*(B3+W*B5)
  387.         fmul.x          (SAVEU,a1),fp0          ; FP0 IS U*V
  388.  
  389.         fadd.x          fp2,fp1                 ; B1+W*(B3+W*B5) + V*(B2+W*B4), FP2 RELEASED
  390.         fmovem.x        (sp)+,fp2/fp3           ; FP2-3 RESTORED
  391.  
  392.         fmul.x          fp1,fp0                 ; U*V*( [B1+W*(B3+W*B5)] + [V*(B2+W*B4)] )
  393.  
  394.         fadd.x          (SAVEU,a1),fp0
  395.         unlk            a1
  396.         rts
  397.  
  398. ;--REGISTERS SAVED FPCR. LOG(-VE) IS INVALID
  399. .pzero
  400.         fmove.s         #0,fp0
  401. .LOGNEG
  402.         unlk            a1
  403.         rts
  404.  
  405. *************************************************************************
  406. * log10():  computes the base-10 logarithm of a normalized input        *
  407. * log2():   computes the base-2 logarithm of a normalized input         *
  408. *                                                                       *
  409. * INPUT *************************************************************** *
  410. *       fp0 = extended precision input                                  *
  411. *                                                                       *
  412. * OUTPUT ************************************************************** *
  413. *       fp0 = log10(X) or log2(X)                                       *
  414. *                                                                       *
  415. * ACCURACY and MONOTONICITY ******************************************* *
  416. *       The returned result is within 1.7 ulps in 64 significant bit,   *
  417. *       i.e. within 0.5003 ulp to 53 bits if the result is subsequently *
  418. *       rounded to double precision. The result is provably monotonic   *
  419. *       in double precision.                                            *
  420. *                                                                       *
  421. * ALGORITHM *********************************************************** *
  422. *                                                                       *
  423. *       log10:                                                          *
  424. *                                                                       *
  425. *       Step 0. If X < 0, create a NaN and raise the invalid operation  *
  426. *               flag. Otherwise, save FPCR in D1; set FpCR to default.  *
  427. *       Notes:  Default means round-to-nearest mode, no floating-point  *
  428. *               traps, and precision control = double extended.         *
  429. *                                                                       *
  430. *       Step 1. Call sLogN to obtain Y = log(X), the natural log of X.  *
  431. *                                                                       *
  432. *       Step 2.   Compute log_10(X) = log(X) * (1/log(10)).             *
  433. *            2.1  Restore the user FPCR                                 *
  434. *            2.2  Return ans := Y * INV_L10.                            *
  435. *                                                                       *
  436. *       log2:                                                           *
  437. *                                                                       *
  438. *       Step 0. If X < 0, create a NaN and raise the invalid operation  *
  439. *               flag. Otherwise, save FPCR in D1; set FpCR to default.  *
  440. *       Notes:  Default means round-to-nearest mode, no floating-point  *
  441. *               traps, and precision control = double extended.         *
  442. *                                                                       *
  443. *       Step 1. If X is not an integer power of two, i.e., X != 2^k,    *
  444. *               go to Step 3.                                           *
  445. *                                                                       *
  446. *       Step 2.   Return k.                                             *
  447. *            2.1  Get integer k, X = 2^k.                               *
  448. *            2.2  Restore the user FPCR.                                *
  449. *            2.3  Return ans := convert-to-double-extended(k).          *
  450. *                                                                       *
  451. *       Step 3. Call sLogN to obtain Y = log(X), the natural log of X.  *
  452. *                                                                       *
  453. *       Step 4.   Compute log_2(X) = log(X) * (1/log(2)).               *
  454. *            4.1  Restore the user FPCR                                 *
  455. *            4.2  Return ans := Y * INV_L2.                             *
  456. *                                                                       *
  457. *************************************************************************
  458.  
  459.         cnop            0,8
  460.  
  461. INV_L10 dc.l            $3FFD0000,$DE5BD8A9,$37287195,$00000000
  462.  
  463. INV_L2  dc.l            $3FFF0000,$B8AA3B29,$5C17F0BC,$00000000
  464.  
  465. *--entry point for log10(X), X is normalized
  466. _log10
  467.         fmove.d         (4,sp),fp0
  468. @log10
  469.         fmove.s         fp0,-(sp)
  470.  
  471.         tst.l           (sp)
  472.         blt.s           .invalid
  473.  
  474.         fcmp.s          #1,fp0
  475.         fbeq            .zero
  476.  
  477.         addq.l          #4,sp
  478.  
  479.         bsr.w           @log                    ; log(X), X normal.
  480.  
  481.         fmul.x          (INV_L10,pc),fp0
  482.         rts
  483. .zero   fmove.s         #0,fp0
  484. .invalid
  485.         addq.l          #4,sp
  486.         rts
  487.  
  488. ;--entry point for Log2(X), X is normalized
  489. _log2
  490.         fmove.d         (4,sp),fp0
  491. @log2
  492.         fmove.x         fp0,-(sp)
  493.  
  494.         tst.l           (sp)
  495.         blt.s           .invalid
  496.  
  497.         tst.l           (8,sp)
  498.         bne             .continue
  499.  
  500.         move.l          (4,sp),d1
  501.         and.l           #$7FFFFFFF,d1
  502.         bne.b           .continue                ; X is not 2^k
  503.  
  504. ;--X = 2^k.
  505.         move.w          (sp),d1
  506.         and.l           #$00007FFF,d1
  507.         lea             (12,sp),sp
  508.         sub.l           #$3FFF,d1
  509.         beq.l           .zero
  510.         fmove.l         d1,fp0
  511.         rts
  512. .zero
  513.         fmove.s         #0,fp0
  514.         rts
  515.  
  516. .continue:
  517.         bsr             @log                    ; log(X), X normal.
  518.         fmul.x          (INV_L2,pc),fp0
  519. .invalid
  520.         lea             (12,sp),sp
  521.         rts
  522.